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<N ; ABSTRACT 
CD , 

^0 ! We consider the growth of a protoplanetary embryo embedded in a planetesimal disk, 

i We take into account the dynamical evolution of the disk caused by (1) planetesimal- 

planetesimal interactions, which increase random motions and smooth gradients in the 
J> ' disk, and (2) gravitational scattering of planetesimals by the embryo, which tends to 

heat up the disk locally and repels planetesimals away. The embryo's growth is self- 
. consistently coupled to the planetesimal disk dynamics. We demonstrate that details 

. of the evolution depend on only two dimensionless parameters incorporating all the 

^SJ , physical characteristics of the problem: the ratio of the physical radius to the Hill 

radius of any solid body in the disk (which is usually a small number) and the number 
of planetesimals inside the annulus of the disk with width equal to the planetesimal 
Hill radius (which is usually large). We derive simple scaling laws describing embryo- 



o 



2 ' disk evolution for different sets of these parameters. The results of exploration in 

c/3 . the framework of our model of several situations typical for protosolar nebula can be 

summarized as follows: initially, the planetesimal disk dynamics is not affected by the 



X 



presence of the embryo and the growth of the embryo's mass proceeds very rapidly in 
the runaway regime. Later on, when the embryo starts being dynamically important, its 
I accretion slows down similar to the "oligarchic" growth picture. The scenario of orderly 

growth suggested by Safronov (1972) is never realized in our calculations; scenario of 
runaway growth suggested by Wetherill & Stewart (1989) is only realized for a limited 
range in mass. Slow character of the planetesimal accretion on the oligarchic stage of the 
embryo's accumulation leads to a considerable increase of the protoplanetary formation 
timescale compared to that following from a simple runaway accretion picture valid in 
the homogeneous planetesimal disks. 

Subject headings: planets and satellites: general — solar system: formation — (stars:) 
planetary systems 



1. Introduction. 



It is now widely believed (e.g. Safronov 1991, Ruden 1999) that terrestrial planets (planets 
like Earth and Venus) have formed by agglomeration of numerous rocky or icy bodies called plan- 
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etesimals (similar to the present-day asteroids, comets, and Kuiper Belt objects). Understanding 
this accretion process is central to understanding planet formation. 

Since the pioneering work of Safronov (1972) it has been known that the dynamics of the 
planetesimal disk, namely the evolution of planetesimal eccentricities and inclinations, is a criti- 
cal factor in this process because relative planetesimal velocities determine to a large extent the 
accretion rate of protoplanetary bodies. The dynamical state of the disk is affected by several 
processes. Gravitational scattering between planetesimals leads to the growth of their cpicyclic en- 
ergy, thereby increasing relative velocities in the disk. This interaction can proceed in two different 
regimes depending on the relation between the random velocity of planetesimals v relative to local 
circular motion and the differential shear across the Hill radius of two interacting bodies rn- 



where mi and m2 are the masses of the interacting planetesimals, Mc is the mass of the central 
star and a is the distance from the central star at which the interaction occurs. Scattering is said 
to be in the shear-dominated regime when v < il,rH, and in the dispersion-dominated regime when 
V > flrn (Stewart Sz Ida 2000), where $7 = ^/GMc/a'^ is the rotation frequency of a Keplerian 
disk. Scattering in the first regime is strong for planetesimals separated by ~ rn and they can 
increase their random velocities very rapidly as a result of it. In the second regime scattering is 
mostly weak even at close encounters because relative velocities are large; in this regime it requires 
much more time to change the kinematic properties of planetesimals. Since scattering in the shear- 
dominated case is so efficient it is likely that this scattering would quickly heat up the disk and 
bring it into the dispersion-dominated regime. For example, a 50-km planetesimal (corresponding 
to a mass ^ 10^^ g) at 1 AU would enter the dispersion-dominated regime for interaction with 
other planetesimals when its epicyclic velocity is as small as 3 m s— 1 (corresponding eccentricity 
is only 10"*^). Thus, it is plausible that at the evolutionary stage when the first massive bodies 
start appearing in the system, planetesimals interact with each other in the dispersion-dominated 
regime — an assumption which we will be using throughout this paper. 

Later in the course of the nebular evolution, another dynamical excitation mechanism emerges: 
gravitational interaction with the newly born massive bodies. We call these bodies embryos (fol- 
lowing Safronov 1972) since they are the precursors of the protoplanets; we will require that their 
masses be greater than the planetesimal masses. When these embryos become massive enough they 
excite planetesimal velocities in their vicinity and also change the spatial distribution of planetesi- 
mals around them. The importance of these effects was first emphasized by Ida &; Makino (1993), 
who used N-body simulations to study planetesimal disks with embedded embryos. 

Depending on the dynamical state of the planetesimal disk, protoplanetary growth can proceed 
in several difi"crcnt regimes. Safronov (1972) argued that emerging embryos heat the planetesimal 
disk up to the point where the random velocities of constituent planetesimals are of the order of the 
escape velocity from the most massive embryos. It is not clear, however, whether the embryo has 
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enough time to increase planetesimal velocities to such large values given inefficiency of scattering 
in the dispersion-dominated regime. Stewart & Wetherill (1988) have noted the importance of 
the dynamical friction in the redistribution of random energy among planetesimals of different 
masses and between planetesimals and embryos. They argued in favor of much smaller planetesimal 
epicyclic velocities which facilitates planetary accretion. However the spatially inhomogeneous 
character of the dynamical effects of the embryos on the disk were overlooked in their study. 

To illustrate the difference between these two scenarios let us estimate the embryo's accretion 
rate dM^/dt in the two-body approximation^ (neglecting the gravity of the central star): 

~ TTRlnmN- 1 + — ^ , (2) 



dt ^ Vz\ ' ReV^ / 

where m is the average planetesimal mass, Mg is the mass of the planet, Re is its radius, N is the 
surface number density of planetesimals, and Vz is the average planetesimal velocity in the vertical 
direction (determining the disk thickness and, thus, the local volume density of planetesimals; 
it is usually of the same magnitude as v). The second term in brackets describes gravitational 
focussing, which can strongly increase the collision cross-section over its geometric value ttR^ when 
the planetesimal velocity is smaller than the escape speed from the embryo's surface ^/2GMe/Re- 
Safronov's (1972) assumption means that focussing is weak. Then one can easily show that 

1 dMe i/'H 

oc M-^/\ (3) 



Me dt 



i.e. the embryo's growth rate slows down as its mass increases. This type of planetary growth is 
known as orderly growth, because it implies that many embryos grow at roughly the same rate. 
On the contrary, if one neglects the embryo's effect on v and and assumes focussing to be strong 
(following Wetherill &: Stewart 1988) one finds that 

1 dMe I/O 

oc M^^, (4) 



Me dt 



i.e. embryo's growth accelerates as its mass increases. This corresponds to the so-called runaway 
accretion regime (Wetherill & Stewart 1989), which allows massive bodies grow very quickly^ in 
contrast with the orderly growth picture. Note that the presence of the dynamical friction is not 
crucial for the onset of runaway growth (it only speeds it up) — what is important is the regulation 
of planetesimal velocities by planetesimal-planetesimal scattering, rather than embryo-planetesimal 
scattering which makes v and Vz independent of Mg. 

However, it is also possible that effects of the embryo on the velocity dispersion are important 
but not so strong as orderly scenario assumes; in this case planetesimal velocities can be increased 



^See also Kokubo & Ida (1996, 1998) for a similar discussion. 

^Formally the condition (4) means that embryo reaches infinite mass in finite time given unlimited supply of 
material for accretion. 
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up to several ^IRh, where Rh = a{Me/Mc)^/^ is the HiU radius of the embryo (note that it is 
considerably larger than the Hill radius rn characterizing planetesimal-planetesimal scattering). 
Then formula (2) yields equation (3) but with a different proportionality constant. This regime is 
called oligarchic (Kokubo &; Ida 1998) and is different from orderly growth because of the role of the 
gravitational focussing (it is very important in this regime) and the values of planetesimal velocities 
(they are much smaller than the embryo's escape speed). As a result, more massive embryos grow 
slower than the less massive ones (similar to orderly growth) but embryos still grow faster than 
planetesimals in their vicinity [similar to runaway growth, see Ida & Makino (1993)]. The accretion 
in the oligarchic regime is considerably slower than in the runaway regime but faster than in the 
orderly regime. For further convenience we briefly summarize major properties of all 3 accretion 
regimes: 

• orderly — growth of velocity dispersion in the disk is dominated by the embryo, gravitational 
focussing is weak; 

• runaway — velocity dispersion growth is determined solely by the planetesimal-planetesimal 
scattering, focussing is strong; 

• oligarchic — velocity dispersion growth is dominated by the embryo, focussing is strong. 

The growth timcscalc of the embryo is a very important quantity. The need for the giant 
planets to accrete their huge gaseous mass onto their rocky cores (core-instability model, see Mizuno 
1980; Stevenson 1982; Bodenheimer &; Pollack 1986) prior to the dissipation of the gaseous nebula 
(HoUenbach et al. 2000) constrains the growth timescale in the giant planet zone — between 5 and 
30 AU in the protosolar nebula — to be < 10^ — 10^ yr. This interval must accommodate several 
distinct evolutionary stages — planetesimal formation from dust, planetesimal coagulation into 
protoplanetary embryos, and finally, catastrophic collisions between embryos to form rocky cores 
of giant planets. Thus, embryos have only part of this time available for their growth. Orderly 
growth is characterized by very long timescales — of order 10^ — 10^ yr (Safronov 1972). On the 
contrary, runaway growth allows embryos to reach ~ 10^^ g (approximately the mass of the Moon) 
in about 10^ yr (Wetherill & Stewart 1993). Runaway growth stops at approximately this mass 
because the embryo accretes all the solid material within its "feeding zone" — an annulus with the 
width of several Rh (Lissauer 1987). This mass is sometimes called the isolation mass Mis (see 
§3.1). 

The dynamical effects of the embryo on the disk, as described numerically by Ida & Makino 
(1993) and Tanaka & Ida (1997) and analytically by Rafikov (2001, 2002a, 2002b, hereafter Papers 
I, II, &; III correspondingly) and Ohtsuki &; Tanaka (2002) are likely to change this conclusion. 
The results of these studies indicate that even if the mass of the embryo is much smaller than the 
isolation mass it will dominate the dynamical evolution of the nearby region of the disk, which 
brings the accretion into the much slower oligarchic mode. In this case the timescale required to 
achieve Mis becomes longer than that predicted by a simple picture of the runaway growth. 
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In this paper we study the coupled evolution of the planetesimal disk and protoplanetary 
embryo. Our goal is to determine (1) which conditions need to be fulfilled for the aforementioned 
regimes to be realized in protoplanetary nebula, (2) whether there are transitions between them 
and when do they occur, (3) which regime sets the timescale for the planetary growth in the disk. 
To answer these questions we consider a very simple model of the embryo-planetesimal disk system. 
Wc describe the behavior of the disk by representing it as a planetesimal population with a single 
nonevolving mass m [similar to "monotrophic" model of coagulation of Malyshkin &; Goodman 
(2001)]. This assumption might seem to be a gross oversimplification but even this toy model can 
still give us a lot of insight into the details of embryo-disk dynamics. The most important omission 
for the dynamics of the system which we make by employing this assumption is the absence of 
dynamical friction within the planetesimal disk. But, as we have mentioned above, this cannot 
preclude the runaway growth from happening. Also, when runaway or oligarchic regimes take 
place, the most massive body in the system grows faster than lower mass planctcsimals because the 
growth rate increases with mass in these accretion regimes, see equations (2) and (4). This justifies 
our assumption of nonevolving planetesimal mass. We devote more time to the discussion of the 
accuracy of our approximations in §4. 

The spatially resolved evolution of the kinematic properties of the disk and the growth of 
the embryo's mass are considered simultaneously and self-consistently within the framework of 
our model. This distinguishes our approach from both conventional coagulation simulations which 
neglect spatial nonuniformitics of the disk properties caused by the embryo's presence (Wetherill 
& Stewart 1993; Kenyon &; Luu 1998; Inaba et al. 2001) and purely dynamical estimates which 
do not allow the embryo's mass to vary^ (Ida &; Makino 1993; Paper I; Ohtsuki Sz Tanaka 2002). 
To describe the coupling between the disk dynamics and embryo's mass growth we use a set of 
equations describing both the planetesimal-planetesimal and the embryo-planetesimal gravitational 
interactions which were derived in Papers II & III correspondingly. Their validity has been checked 
elsewhere using direct integrations of planetesimal orbits in the vicinity of the embryo (Paper III). 
These analytical equations are solved numerically in §3. In addition, to highlight the importance 
and better understand the behavior of different physical processes operating in the system we 
provide simple scaling estimates for the evolution of the disk kinematic properties and the growth 
of the embryo's mass in §2. We discuss the applications and limitations of our analysis in §4. 

We will often use the model of the Minimum Mass Solar Nebula (MMSN, Hayashi 1981) to 
obtain typical numerical values of physical quantities. In particular, we will use the following 
distribution of the surface mass density of the solids in MMSN (Hayashi 1981): 

S = 20gcm-2 (3^)"'^', (5) 
which is obtained assuming that solids constitute 1% of the nebular mass. 



^Probably the closest analog of our method is the multi-zone coagulation simulations of Spaute et al. (1991) and 
Weidenschilling et al. (1997). 
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2. Scaling arguments. 

In this section we provide simple scaling-type solutions describing the coupled evolution of 

the embryo-planetesimal disk system. The reader interested only in the numerical results of more 
detailed and accurate calculations can skip this section and proceed directly to §3. 

We first derive a set of evolution equations for the embryo's mass growth and for the planetes- 
imal velocity dispersions in the horizontal and vertical directions using qualitative reasoning. This 
will allow us to better understand different dynamical mechanisms operating in the disk and predict 
the dependence of the system's behavior on its basic parameters. The equations are derived for 
two different regimes of the embryo-planetesimal gravitational interaction — dispersion-dominated 
and shear-dominated, although we mostly concentrate our attention on the former. We will assume 
that planetesimal-planctcsimal interactions always occur in the dispersion-dominated regime which 
almost certainly should be true. From now on by different velocity regimes we will imply different 
regimes of the embryo-planetesimal interaction. 



2.1. Derivation of the qualitative equations. 

We start with the shear-dominated regime. In this case the approach velocity of a planetesimal 
towards the embryo is ~ URh and the effect of the tidal force of the central star is quite important. 
As a result scattering is strong and each encounter with the embryo imparts the planetesimal with 
a change of its velocity in the horizontal direction Av ~ CIRh (Petit & Henon 1986). At the 
same time the excitation of the vertical component of the planetesimal velocity is weak, because 
the thickness of the planetesimal disk is smaller than the embryo's Hill radius Rh- As a result, 
scattering is almost two-dimensional and scattered planetesimals are restricted to the plane of the 
disk. Only the deviations from purely planar geometry can lead to the growth of inclination. On 
the basis of geometric considerations we can predict that Avz ~ {vz/^Rh)^v ~ Vz- 

Planetesimals at different semimajor axes separations have different synodic periods (the time 
between successive conjunctions with the embryo). In the shear-dominated regime the most impor- 
tant region of the disk is that within a Hill sphere of the embryo. Planetesimals in this region have 
radial separations from the embryo ~ Rjj and travel past the embryo at relative velocity ~ Q.RH- 
Thus, the typical time between consecutive conjunctions with the embryo for planetesimals in this 
part of the disk is 



r_ = ljs^-i-^~j^-i-^. (6) 

^ Z Rh Rh 



Then we can approximately describe the dynamical evolution of the planetesimal disk due to the 
embryo in the shear-dominated regime as 

n^iQRnf, (7) 



dt 



sd 
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dt 



em— pi 



(8) 

film. 



We now turn to the dispersion-dominated regime, first focussing on embryo-planetesimal scat- 
tering. Here the two-body approximation usually proves to be very useful. In this approximation 
the influence of the third body — the central star — can be neglected and three-body scattering 
problem naturally reduces to the two-body problem which is tractable analytically. Wc will carry 
out our qualitative arguments using the instantaneous scattering approximation. In this approxi- 
mation a body with a mass mi approaching another body with a mass m2 with relative velocity v at 
an impact parameter b experiences a change in velocity roughly equal to Av ~ Gm2/{bv) assuming 
small angle scattering. This approximation is valid if 6 > 60 ~ G{m\ + m-2)lv^^ since for 6 < 60 
large-angle scattering takes place. The upper bound on b is set by the smaller of the disk vertical 
thickness H ~ v^j^ ^ ^0 ™d the range of planetesimal epicyclic excursions vjVl (usually H is 
the one). Since in three-dimensional scattering every decade in the impact parameter h between 
60 and B. provides roughly the same contribution to the scattering coefficients we find, integrating 
the expression for Au over 6, that the change of per encounter averaged over b is 



(Av^) K In A, where A «i H/bo ~ 77^ — • (9) 

\ V J vvz (ila)'^ m\ + m2 

Changes of and v1 are different only by constant factors because scattering has a three-dimensional 
character in the dispersion-dominated regime. 

Embryo-planetesimal scattering leads to changes of planetesimal velocity only at conjunctions 
with the embryo. Planctcsimals which can experience close encounters with it occupy an annulus 
of the disk with width ~ around the embryo's orbit. The typical time between the consecutive 
approaches to the embryo for the planetesimals in this zone is ~ Tsyni^Rn/v)- Then we can write 
the following equation for the velocity evolution caused by the embryo-planetesimal encounters in 
the dispersion-dominated regime: 



Tsyn ^Rh a V^Vz 

da 



where Ag is A computed using m2 = Mg S> m = mi. 

We next turn to planetesimal-planetesimal scattering, which always occurs in the dispersion- 
dominated regime (see §1). The volume number density of planetesimals n can be expressed 
through the surface number density as n ~ N/H ~ NQ/vz- Then we can write the rate of 
random velocity growth in the disk due to planetesimal-planetesimal encounters as 



d / 2 2\ 



dt 



pi~pi 



„2 i^rn)^ 



vn-^{Av') « ^Nrjj '^ InA^, (11) 

i VVz 
dd 
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where m is the typical planetesimal mass and Ap is the value of A calculated using mi = m2 = m. 
More accurate forms of the velocity growth equation (11) in the dispersion-dominated regime have 
been previously obtained by a number of authors (Stewart k, Wetherill 1988; Ida 1990; Stewart 
&; Ida 2000; Paper II). One can compare them with (10) and see that our equation really is a 
qualitative version of more accurate results. 

The embryo's accretion rate can also be approximated by the two-body expression (2). We 
will assume that the physical collisions between the embryo and planetesimals occur with velocities 
at infinity considerably smaller than the embryo's escape velocity. This assumption, which implies 
that gravitational focussing plays an important role in the accretion process, is usually justified in 
the early stages of the embryo's growth but can run into problems later if planetesimal disk heats up 
very strongly. We will later check the validity of this assumption (see §4). In the shear-dominated 
regime the approach velocities of planetesimals are dominated by the differential shear in the disk, 
i.e. V ~ Q,Rh- In the dispersion-dominated regime the approach velocities are dominated by the 
epicyclic motion of planetesimals. As a result we find that 



M 

in the shear-dominated regime, and 

M 



^ mN^ (12) 

sd Vz 



^„,N^^R^ (13) 

dd VVz 



in the dispersion-dominated regime. More accurate expressions for the accretion rate in these 
regimes can be found in Greenzweig &; Lissauer (1992) and Dones & Tremaine (1993). 

At this point it is useful to switch to the following dimensionless notation 

,2_ jv') ,2_ {vD j^^M^ R^ 
{nrnr ' i^rny' rn ' ^ R„' 

T = t/tsyn, where tgyn = ^O"^— = Tgyn— = TsynM^f^, (14) 

6 rn rn 

i.e. we normalize the planetesimal velocities by the shear across the planetesimal Hill radius rn = 
a{m/Mc)^/^ , time by synodic period at 1 rn and embryo's mass by the planetesimal mass m. The 
parameter p is the ratio of the embryo's physical radius to its Hill radius, which remains constant as 
the embryo grows at constant density. Note that in this notation embryo-planetesimal interaction 
is in the shear-dominated regime when s < M^^^, and in the dispersion-dominated regime when 

S > MV3. 

With these definitions we can rewrite the equations governing the behavior of the planetesimal 
velocity and embryo's mass growth in the following dimensionless form: in the shear-dominated 
regime (s < M^^^) 
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^^iNarH)^^ + slM'/', (16) 

CtT z 

-^^p{NarH) — . (17) 

dT Sz 

In the dispersion-dominated regime (s > A^^/^) 

^ {s\ .^) ^ (iVar^)^ + Al^^, (18) 
dr ^ ' ssz s^'Sz 

dM M'^/^ 

_»p(iV„.„)--. (19) 

The different factors entering the Coulomb logarithms are 

Ap = s^Sz and Ag = s^Sz/M. (20) 

In the equations of evolution of ,s and (15), (16), and (18) the first term on the r.h.s. de- 
scribes the effect of planetesimal-planetesimal scattering while the last one accounts for the embryo- 
planetesimal encounters. We have dropped all numerical factors entering these equations because 
our qualitative analysis cannot determine them. However this should not affect our general con- 
clusions. 

We do not write down an equation of evolution for N — the surface density of planetesimals. 
Instead we simply assume it to be constant. This is certainly justified while the embryo's mass is 
small and it cannot perturb the spatial distribution of planetesimal orbits around it. However as it 
grows bigger, the embryo starts repelling planetesimals (Papers I &: III) which affects their surface 
density. This effect is not included in our approach. It is often argued that the conservation of 
the Jacobi constant of planetesimals in the course of their scattering by the embryo can maintain 
their instantaneous surface density (which matters for the accretion) at more or less constant level; 
unfortunately it will turn out (see §3) that this assumption is not very good for very large embryo 
masses and more accurate treatment of evolution of N is required. However, this disadvantage of 
the assumption N = const does not affect our major conclusions. 

More accurate analysis shows that as well as the heating terms the equations for the plan- 
etesimal velocity evolution should also contain transport terms. These appear because mutual 
gravitational scattering of planetesimals acts as a source of effective viscosity in the disk (Papers I 
& II; Ohtsuki & Tanaka 2002) which tends to smooth out any inhomogeneities of planetesimal sur- 
face density and velocity dispersions. Equations (15)-(18) do not take this effect into account but 
this does not degrade the accuracy of our conclusions: it follows from our analysis of planetesimal- 
planetesimal scattering in Paper II that the transport (or viscous) terms lead to the effects of the 
same magnitude as the planetesimal heating terms. We will discuss this issue in more detail in §4. 

The system (15)-(20) describes the dynamics of the disk only in the immediate vicinity of 
the embryo — within ~ Rh from the embryo's orbit in the shear-dominated regime and within 
~ sRh / M^^^ in the dispersion-dominated regime. We will call this region the ''heated" zone. Only 
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the dynamics of this region matters for the embryo's growth because embryo can accrete material 
only from this part of the disk. Outside the heated region, the embryo's influence is negligible and 
disk evolution proceeds almost exclusively under the action of planetesimal-planetesimal scattering. 
Evolution of the distant parts of the disk is described by equation (18) with the last term in its 
r.h.s. dropped and is rather simple. 

One can easily see that the evolution equations incorporate only two free parameters char- 
acterizing disk properties: Naru and p. What is more important, the more accurate equations 
derived previously in Papers II &: III also exhibit this property (see §3 and Appendix A). This 
implies that various evolutionary scenarios describing the growth of the embryo can be classified 
completely by these two parameters. After such classification is performed the general features of 
the evolution can be predicted once the values of the parameters Narn and p are known. Typical 
values of these parameters in MMSN are discussed later in §3. 

Our subsequent analysis of equations (15)-(19) will be essentially asymptotic: we will separate 
the embryo-disk evolution into a sequence of regimes in which only some terms in the evolution 
equations are important. Evolution in the transition zones between separate asymptotic regimes 
will be neglected. The unimportance of the transient stages typically requires the changes of various 
quantities in asymptotic regimes to be much larger than in the transition zones. We will see that 
this is usually true only when extreme values of parameters Narn and p are chosen. However, even 
for a realistic choice of these parameters typical for MMSN our analysis will still give us important 
insight into the coupled evolution of the embryo-disk system. 



2.2. Separation of the velocity evolution regimes. 

Planetesimal-planetesimal scattering dominates the evolution of the planetesimal velocity when 
the first term in the r.h.s. of the corresponding equation exceeds the second term representing the 
effect of the cmbryo-planetesimal scattering. The shear-dominated regime corresponds to > 
in s — M coordinates. One can see from (15) that in this regime mutual planetesimal encounters 
dominate the growth of s and Sz when 

M<Msd{s)^\nAp^^^. (21) 

In writing down this condition we have assumed that s ~ which is always the case for three- 
dimensional scattering characterizing planetesimal-planetesimal encounters. Note that when s ~ 1, 
i.e. when even planetesimal-planetesimal scattering is close to the shear-dominated regime, the 
mass at which the embryo begins to dominate the disk dynamically is Mg ~ m{NarH) [InAp ~ 1 

when s,,Sz ~ 1, sec (20)]. The same estimate of Mg was obtained in Paper I where both mutual 
planetesimal and embryo-planetesimal encounters were treated in the shear-dominated regime. 

The dispersion-dominated regime occurs when M < s^. In this case the disk dominates its 
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Fig. 1. — Different regions of embryo-planetesimal interactions, in coordinates s = v/{Q.rH) and 
M. = Me/m [see equation (14)]. The thin sohd hne separates shear- and dispersion-dominated 
regimes, while the dashed curve separates regions where either embryo or planetesimals dominate 
disk heating [see equations (21) and (22); deep in the dispersion-dominated regime it has a slope of 
1/2]. The embryo controls the disk evolution in the shaded region, mutual planetesimal scattering 
dominates disk evolution in the unshaded region. Two typical evolutionary tracks for an embryo 
of fixed mass are displayed: for M.a < -M* and for > -M.*- 

own dynamical evolution only when [see (18)] 

M < Maais) - (jj^) (iVar^)V2.V2. (22) 

A similar condition (but without InAg) was previously derived by Ohtsuki &; Tanaka (2002). 

Conditions (21) and (22) separate the s-M plane into regions where disk self-heating dominates 
over the effect of the embryo, and vice versa. These regions are displayed in Figure 1 and the one 
where embryo controls the disk heating is shaded. The curves in the s-M plane given by (21) and 
(22) intersect with each other and with = at the point (s*, A1*) where 



^ {InApf/^NarHf/^ ^ (In Aj,)V5(iVari/)^/^ = mI^^ Ap ^ {Narnf^ (23) 
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where we have taken into account that In Ag ~ 1 when s, Sz ~ Ai^^^, corresponding to the transition 
between the shear- and dispersion-dominated regimes [see (20)]. 

Let us assume initially that the embryo's mass docs not evolve and the disk is initially in the 
shear-dominated regime. Since the disk can only heat up, the evolutionary track of the embryo- 
disk system in the s-M. plane is a straight line going to the right at constant M. It then follows 
from Figure 1 that when M < Mi,, the embryo can never be effective at heating the disk up 
and the evolution of eccentricities and inclinations of planetesimals takes place only as a result of 
planetcsimal-planctcsimal gravitational scattering — the effect of the embryo-planetesimal inter- 
action is much weaker (evolutionary track A- A' in Figmc 1). On the contrary, when Ai > 
(track B-B' in Figure 1) embryo will start dominating the dynamical evolution of the heated zone 
when the disk is sufficiently excited in the shear-dominated regime (at point b) . As the heated zone 
gets more and more excited (s grows) mutual planetesimal scattering starts eventually to dominate 
the disk evolution again (at point b'). This is actually not surprising because one can see from 
(18) that in the dispersion-dominated regime the embryo's contribution to the disk evolution scales 
roughly as while that of the mutual planetesimal encounters behaves as at a fixed Ai. 
Thus, as s increases planetesimal encounters finally take over the control over the disk evolution 
in the heated region from the embryo. It would however take a very long time for the system to 
reach the transition point 6' because usually Sb> 3> 1 and disk random motion evolution is very slow 
in the dispersion-dominated regime. We will not discuss any other details of the disk dynamical 
evolution with fixed M because realistic embryos have p / and their masses have to increase. 



2.3. Separation of the mass and velocity evolution regimes. 

Whenever the mass of the embryo is allowed to vary one has to study the coupled embryo- 
disk evolution using the full set of equations (15)-(17) or (18)-(19). We will concentrate here on 
exploring the dispersion-dominated regime (we briefly touch on the shear-dominated regime in the 
end of this section). 

Let us estimate different timescales associated with the evolution in the dispersion-dominated 
regime. Since Sz ~ s, the timescale of disk evolution caused by planetesimal-planetesimal scattering 
is [first term on the r.h.s. of (18)] 



Tpi 



1 ds 
s dr 



i 4 1 

^ S 1 

• 

pi ) In Ap NavH 



(24) 



The timescale characterizing the effects of the embryo is given by [last term on the r.h.s. of (18)] 



s 



Ids ^ "5 



s dr 



In Ae M. 



(25) 



Finally the characteristic time of the embryo's mass growth is 



1 dM\-'^ .2 1 
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Fig. 2. — Separation of different regions in the s-Ai space [see equation (14)] in the case of coupled 
evolution of disk velocity dispersion and embryo's mass. The dotted line displays the asymptotic 
relation (28) and the short-dashed line shows relation (27). The dot-dashed line represents the 
restriction put by the existence of the isolation mass [see (38) in §3.1]. The embryo's mass evolves 
faster than the disk velocity dispersion in the shaded region (and vice versa in the unshaded region). 
The thick solid curve shows typical evolutionary track of the embryo-disk system. The long-dashed 
curve has the same meaning as in Figure 1 (its slope is 1/2 when s 3> A4^^^). 



By requiring Tpi = Tg^ we retrieve the separation condition (22) again. 

Let us assume that at the beginning embryo has mass M < Mi, [Mi, is given by equation 
(23)]. Then initially disk evolution is dominated by planetesimal-planetesimal scattering. The disk 
evolves faster than Me grows only when Tpi < Tm or when (see Figure 2) 

M<MM-A^)^(^^y ^. (27) 

When M grows bigger than Mm given by (22), the embryo starts dominating the velocity evolution 
of the disk in the heated zone. This evolution proceeds faster than the embryo's mass grows when 
Tem < Tm or 

3/5 

s^/\ (28) 



M > MM-em{s) 



p{NarH) 
InAe 



- 14 - 



The curves in the s-Ai plane given by equations (27) and (28) intersect with each other and with 
(22) in the dispersion-dominated regime (when s > Ai^^^) if p > po, where 

Thus, when p > pQ there is a region in the dispersion-dominated regime bounded by curves 
(27) and (28) in which the embryo's mass grows faster than the random energy in the heated zone 
(shaded region in Figure 2). Everywhere outside this region the disk heats up faster than the 
embryo grows. This allows us to sketch the evolution of the system in the following way (thick 
solid curve a — e in Figure 2) starting with some Mo < Mm -pi- initially the disk heats up faster 
than Ai grows which produces an almost horizontal track a — 6 in the s-Ai plane — see Figure 2. 
As s gets to the point b where Aio ~ Mm-pU mass starts growing faster than s; as a result the 
evolution track becomes almost vertical {b — d). At point c when the condition Ai = Mdd [see (22)] 
is fulfilled, the embryo takes control over the dynamical evolution of the heated zone but its mass 
still grows faster than s increases [Tm < Tem-, see (25) and (26)]. This continues until evolutionary 
track reaches the curve M = Mm -em at point d. At this point timescales of the embryo's growth 
and disk dynamical evolution become almost equal and embryo-disk system evolves stably along 
the curve d — e. 

This general picture is confirmed by the solution of the system (18)-(19). When A4 < Mdd a-nd 
planetesimal-planetesimal encounters dominate the disk heating one finds that (treating logarithmic 
factors as constants) 

s{t) « [{NavH) In Ap]^/VV4, M{t) ^Mo[i- (r/r^^)^/'] , (30) 

where 

Ml^\NarH)p^ 

is the timescalc of the runaway growth. One can see from this solution that when the embryo is not 
massive enough to affect the disk dynamics around it does grow in the runaway regime which is in 
complete agreement with the scenario of Wetherill & Stewart (1989, 1993). However, in practice 
the runaway solution (30) is only valid until A^(r) « Mdd- This mass is reached when r Tmn 
since typically Mdd S> Mq. This means that embryo starts to dominate disk heating at (see Figure 
2) 

1/2 



Si = s{Trun) ~ ( ^" ) , Ml = Mdd{si) 

\pMrJ 



(iVari^) (In Ap)3/2 



p 



V27Wy^ InAe 



(32) 



After this happens evolution proceeds along a different route. One can find from (18)-(19) 
that for 2> All the evolutionary track is described by the equation 



M 



(33) 
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Asymptotically, when s 3> si this relationship reduces to (28) and one obtains that 

^' (lnAe)V7 ^ ' (lnAe)6/7 ^ ■ ^^^^ 

Note that this solution does not depend on M.q, Mi, s\, etc. — all the memory of initial conditions 
is lost. For this reason we would expect any evolutionary track describing embryo-disk system to 
behave eventually as equation (34) predicts, independent of the initial conditions used. 

The solutions of evolution equations obtained in this sections directly pertain only to the case 
p < Pq. If parameter p is such that p > po the curves given by (27) and (28) do not intersect in 
the dispersion-dominated regime. As a result embryo-planetesimal system can in principle switch 
from this velocity regime to the shear-dominated or the intermediate velocity regime. We do not 
explicitly study the case p > po in this paper — it turns out to be not very different from the 
case p < Po for the typical MMSN parameters (see §3). We also do not derive any results for the 
embryo-disk evolution in the shear-dominated regime. There are several reasons for this: 

• One would expect that initial disk velocity dispersion and starting mass of the embryo are 
correlated with each other since both require some time to grow; this is likely to bring the 
system into the dispersion-dominated rather than shear-dominated regime. 

• For the parameters typical for MMSN the variation of s in the shear-dominated regime is 
restricted to a rather narrow range. It is thus possible that the asymptotic analysis such 
as presented in §2.2 and §2.3 would not be very helpful in this case since the effects of 
the transition zones between the different asymptotic regimes will be quite important (see 
discussion in the end of §2.1). This is even more important for the case p > pq. 

• It is unlikely that the embryo-disk system can spend a long time in the shear-dominated 
regime since the dynamical evolution of the disk is usually very rapid in this case. Thus the 
shear-dominated evolution is expected to be only a transient stage of the system's history. 

• The numerical approach we are going to use for the description of the embryo-planetesimal 
interaction is not very accurate in the shear-dominated regime (see Paper III). 

These arguments provide the basis for our present neglect of the shear-dominated and inter- 
mediate velocity stages of the embryo-disk evolution. However it is clearly not difficult to study 
these regimes qualitatively based on equations (15)-(17) in the spirit of §2.2 and §2.3 if the need 
arises. 



3. Numerical results. 



One can study the embryo-disk interaction and check the simple qualitative predictions of §2 
using more robust machinery. To do this we use the apparatus elaborated in Papers II &: III for 
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the statistical description of the planetesimal-planetesimal and embryo-planetesimal gravitational 
scattering. There we have derived a set of equations governing self-consistently both the planetesi- 
mal disk dynamics and the embryo's mass growth, taking into account planetesimal-planetesimal as 
well as embryo-planetesimal gravitational scattering. These equations are the most accurate when 
the embryo-planetesimal interaction is in the dispersion-dominated regime although wc arc able to 
provide an approximate treatment of the intermediate velocity regime as well (sec Paper III). In 
Appendix A we transform this set of equations into a dimensionless form. We then numerically 
solve them and find the embryo's mass as a function of time; in addition we obtain time-dependent 
and spatially-resolved profiles of the planetesimal surface density and dispersions of eccentricity 
and inclination. What makes this analysis novel is that it takes into account the effects of plan- 
etesimal disk nonuniformities induced by the embryo's gravity. The closest existing analog of our 
method is a multi-zone coagulation simulation approach of Spaute et al. (1991) and Weidenschilling 
et al. (1997) which achieves spatial resolution by employing a large number of single-zone coagu- 
lation simulations interacting through the boundaries. However in their published simulations the 
treatment of the disk dynamical evolution is very rough and the degree of the spatial resolution is 
not very large — of the order of Hill radii of biggest bodies at best while our approach allows us 
to use much smaller grid sizes. 



3.1. Numerical setup. 

In this paper we are interested in the growth of a single isolated embryo. We assume a 
planetesimal population with a single characteristic mass and disregard the evolution of this mass 
due to the planetesimal coagulation. We neglect any mechanisms which can affect planetesimal 
dynamics (gas drag, inelastic collisions, etc.) other than gravitational scattering. The embryo's 
eccentricity and inclination are assumed to be always zero. We are going to relax these assumptions 
in future work. As the embryo's mass increases, its Hill radius Rh (which is a characteristic 
length scale of the problem) grows as well. For this reason we compute the spatial distributions 
of various quantities not in the physical coordinates but in Hill coordinates of the embryo which 
naturally adjust to the embryo's mass growth. For the same reason in our numerical calculations 
we characterize the dynamical state of planetesimal disk not by s and Sz defined by (14) but by S 
and Sz — velocity dispersions in horizontal and vertical directions normalized by the embryo's Hill 
radius Rh (and not by the planetesimal Hill radius rn)' 

{^Rh? ^2/3 _y^2/3' {^Rnf f/J^ M2/3' ^^^^ 

where CTg and CTj are the (nonreduced) eccentricity and inclination dispersions of Paper II. When 
making comparisons with §2 we usually switch back to variables s and Sz- The dimensionless vari- 
ables A4 and r defined in equation (14) are used instead of the embryo's mass and time. The width 
of the integration region in Hill coordinates is chosen to be large enough [usually (—AORh', 4:0Rh) 
in coordinates which adjust to the embryo's mass] so as to minimize boundary effects. 
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To provide an estimate of the random velocity dispersion s in the heated zone (which is one 
of the functions of interest) in our numerical solutions we have chosen the following method: we 
determine the extent of the region around the embryo within which the excess of eccentricity 
dispersion over its value at infinity is > 20% of the maximum value of this excess. We define this to 
be the heated zone. Then we average (S"^ + 5*^)^/^ over this region and use the resultant quantity 
[multiplied by Ai'^^^, sec (35)] as a measure of planetesimal random velocity s. At this point we do 
not make a distinction between the total velocity ys^^Tsf and its horizontal component s because 
they are approximately the same at our level of accuracy. This is rather loose definition of s but it 
suffices for our purposes. 



As we have already mentioned in §2.1 the outcome of the embryo and planetesimal disk evolu- 
tion depends only on the values of Narn and p. This simplifies our treatment a lot by combining 
several different quantities (such as m, a, N, M^, etc.) into just two parameters. The parameter 

Narn has the meaning of (roughly) the number of planctcsimals within the annulus of width rib- 
and radius equal to the embryo's scmimajor axis a. This number is typically quite large. Parameter 

p is simply the ratio of the physical size of the body to its Hill radius. Since both radii depend 

1/3 

on the mass of the body in the same way (both are oc Mg ) p is actually independent of Me and 
can be taken to characterize both the embryo and planetesimals. In many astrophysically relevant 
situations this parameter is small. Using MMSN parameters given by (5) we can estimate that 

= ^i^JmF" \~) W * ' 



3 McV^^ ^^_3 /I AU\ /3 g cm-3\ /^^^X 1/3 



and 

where p is the density of solid material constituting planctcsimals as well as the embryo. In Table 
1 one can find typical values of these parameters in different important locations of the protosolar 



Table 1. Typical values of Narn and p} 



Typical object 


a, AU 


Narn {m = 10^^ g) 


p 


Earth 


1 


360 


3.6 X 10-3 


Jupiter core 


5 


800 


7 X 10"^ 


KBO 


40 


2.3 X 10^ 


0.9 X 10-^ 



^In computing these values M^, = Mq and p = 3 g cm 
are used. 
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nebula. Note that for a given p and equation (37) unambiguously determines the semimajor 
axis in the MMSN corresponding to a particular value of p. Knowledge of Narn then allows one 
to fix the planetesimal mass m. 

We numerically solve the evolution equations of Appendix A and compare the results with 
the predictions of §2. We display the outcomes of calculations for 3 representative cases: Narn = 
400, p = 0.004 corresponding to a = 0.9 AU and m = 7.3 x 10^° g which is typical for the terrestrial 
planet region; Nar^ = 10^, = 10~^ (a = 3.6 AU and m = 1.5 x 10^^ g) which is close to the 
region of giant planet formation; Narn = 10^, p = 10^^ (a = 36 AU and m = 4.7 x 10^^ g) 
which is representative of the distant part of the protosolar nebula where Kuipcr Belt Objects 
(KBOs) arc assumed to form (although the typical planetesimal mass is almost certainly too big 
in this case). In each case the calculation was carried out for several values of initial embryo's 
mass M.Q and velocity dispersion sq (which is usually unimportant). Each calculation starts with a 
protoplanetary embryo in a homogeneous disk; we follow both the subsequent growth of embryo's 
mass and the development of disk inhomogeneities. 

We allow the disk-embryo system to evolve until M exceeds the restriction put by the existence 
of the isolation mass Mis (for higher masses our results are not accurate because planetesimal surface 
density would be considerably reduced by accretion). From the definition of Mis (see §1) it follows 
that^ Mis ~ 27ra x 2a{Mis/Mc)^/^mN , or in our dimensionless notation 

Mis - (47r)3/2(iVarij)3/2. (38) 

In fact, in the dispersion-dominated regime this condition should be additionally multiplied by 
S^^^ = s^^'^/Ai^^^ because in this velocity regime the width of the feeding zone is about the 
planetesimal epicyclic excursion SRh- Thus, (38) should be considered a lower limit on the isolation 
mass which is valid in the shear-dominated regime only; in the dispersion-dominated regime the 
isolation mass becomes 

Mis - ^■K{NarH)s. (39) 

Obviously the more dynamically excited planetesimal disk is (i.e. the higher s is) the more massive 
the embryo can grow, because it can accrete planetesimals from a larger region of the nebula. 

3.2. General description of evolution. 

The results of the numerical solution of the evolution equations are presented in Figures 3-12. 
We use the following general notation when plotting evolutionary tracks of different quantities {M 
or s): tracks corresponding to the regime where planetesimal-planetesimal scattering dominates 
the disk dynamics are plotted in blue, while those in the regime dominated by embryo-planetesimal 



''We Eissume here that the width of the feeding zone is 2Rh rather than Rh as has been cissumed in Paper I. 
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Fig. 3. — Evolution of the embryo-disk system in s-M coordinates for Narn = 400 and p = 0.004 
(typical for the Earth-forming region of MMSN). The evolutionary tracks consist of two portions: 
where planetesimal scattering dominates the dynamical evolution of the disk (blue dots) and where 
embryo dominates this process (red dots). Curves separating various regions in this plot have the 
same meaning as in Figure 2. The thin solid curves initially coinciding with the dotted tracks 
describe the runaway growth of "passive" embryos. 

scattering are shown with red. The separation between these two regimes is determined based on 
Figures 10-12 (see below). We also plot on each Figure (on the right and top borders) the numerical 
values of physical variables Mg (mass of the embryo) and v (linear velocity of planetesimal epicyclic 
motion) corresponding to their dimensionless analogs M and s. In doing this we take into account 
the aforementioned property [see discussion after (36) and (37) of Narn and p to unambiguously 
specify a and m when p and Mc are fixed. Also, where appropriate we plot not only the evolutionary 
tracks including the effects of embryo's gravity on the disk dynamics ("active" embryo, colored dots) 
but also the curves calculated assuming that embryo does not affect the disk dynamics (thin solid 
lines); in other words we allow the disk to stay completely homogeneous and to heat up only 
due to the planetesimal-planetesimal scattering, while still allowing the embryo to grow ( "passive" 
embryo). This would correspond to the outcome of conventional coagulation simulations (in our 
problem setting) which do not take embryo- induced disk inhomogeneities into account. 




In Figures 3-5 we plot the evolution of the embryo-disk system in coordinates s-M for different 
values of Narn and p. This can be compared to the schematic calculations in Figure 2. Each of the 
plots displays several evolutionary tracks shown by sequences of dots. These tracks usually start 
with moderate values of TWq (from 10 to 10^) and rather small values of sq; so that the embryo-disk 
interaction can start either in the shear- or in the dispersion-dominated regimes. A brief inspection 
of Figures 3-5 shows that the coupled evolution of the embryo-disk system can be naturally split 
into three distinct regimes. In the first regime the mass of the embryo stays almost constant^ and 
only the velocity dispersion in the disk grows. After a short intermediate stage when both s and 
M evolve at comparable rates, the embryo grows faster than disk heats up because of the fast 
runaway accretion. At some intermediate point during this stage the embryo begins to dominate 
the dynamical evolution of the planetesimal disk (blue dots change to red dots). As a result the 
evolutionary track gradually departs from the runaway path and accretion proceeds at a slower 
rate. Finally, the system settles onto the asymptotic regime in which embryo's mass grows as a 
power law of the velocity dispersion in the heated zone s. 



^This happens because the timescale of the runaway growth is longer than the timescale of the disk heating. 




Fig. 5. — The same as Figure 3 but for Narn = 10^ and p = 10 ^ (typical for the KBO-forming 
region of MMSN). 

Figure 6 provides more insight into what is going on in the disk by plotting the time evolution 
of the spatial distributions of planetesimal surface density N and eccentricity and inclination dis- 
persions S and Sz [scaled by Rh, see (35)]. These distributions are plotted for an evolutionary track 
with Narn = 400, p = 0.004 and initial values of mass and velocities A^o = 10,50 = •S'zo = 0.5. 
We also plot the locations of the points at which the snapshots of spatial distributions are taken 
in Figures 3, 7, and 10 (black stars on the corresponding evolutionary track). One can see that 
initially the disk heats up very quickly: S grows from 0.5 to 2.8, Sz from 0.5 to 1.6 on a rather short 
timescale. The ratio of Sz/S is close to 0.55, which is characteristic for homogeneous Keplerian 
disks (Ida et al. 1993; Stewart & Ida 2000). One can see that the embryo's effects are negligible 
and disk is very uniform (red curve). After that rapid runaway accretion of the embryo leads to the 
decrease of S and Sz because growth of Rh surpasses the increase of s and Sz- Disk inhomogeneities 
are rather small while disk dynamics is still dominated by the mutual planetesimal encounters (blue 
curve). However, as soon as the embryo takes over the control of the disk heating (green curve) 
deviations from the uniform state become important: sharp features appear in the distributions 
of S and Sz and the embryo carves out a gap in surface density of planetesimal guiding centers 
around itself (results analogous to those of Papers I &: III). After the intermediate stage (magenta 
curve), the disk finally settles into the asymptotic regime (black curve) and the difference between 
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Fig. 6. — Time evolution of the spatial distributions of planetesimal surface density N (normalized 
by its value at infinity), eccentricity and inclination dispersions S and Sz scaled by the embryo's Hill 
radius for Narn = 400, p = 0.004, AIq = 10, 5*0 = 0.5, Szo = 0.5. Curves of different color represent 
the snapshots of the distributions at different moments of time marked with corresponding color 
(also marked on the evolutionary track in Figures 3, 7, & 10 with starred dots). The dotted line 
shows the initial state of the disk. See text for more details. 



the planetesimal velocities in the heated zone and far from the embryo steadily grows with time. 
The detailed shapes of the spatial distributions of different quantities are almost invariant in the 
asymptotic regime — we would obtain the same final profiles if we were to start with a different 
set of ^AQ^ SeOT SiQ. Note that far away from the embryo, the disk always stays uniform and Sz/S 
is always near 0.55^. 

It is interesting that for a rather long period of time disk evolution in the heated zone takes 
place almost in the shear-dominated regime. This happens because the runaway growth of the 
embryo's mass is much faster than the growth of random velocities in the disk. This switches 
embryo-planetesimal interactions from the dispersion-dominated into the intermediate velocity 



^This ratio of Sz/S eventually falls in the very late stages of embryo-disk evolution, because eccentricity is rather 
effectively excited by the distant encounters (see e.g. Hasegawa & Nakazawa 1990). 
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Fig. 7. — Growth of the embryo's mass in the case Narn = 400, p = 0.004 for different initial 
conditions. Mass as a function of time r is plotted by dots. The thin solid lines display 
corresponding runaway curves obtained by neglecting embryo's dynamical effects. The dimensional 
values of time t and mass of the embryo Me are obtained in the same way as in Figures 3. The 
dot-dashed line represents the isolation mass in the shear-dominated regime defined by (38) and 
by (39) in the dispersion-dominated regime (see text). The dashed line has a slope of unity and 
shows the asymptotic behavior of M. [see equation (45)]. 

regime. Some minor details of the spatial distributions in this regime might be not real because 
our analytical apparatus developed in Paper 111 is not very accurate in the intermediate velocity 
regime; however the overall description of the disk evolution should be robust. 

We have mentioned in §2.3 that planetesimal-planetesimal scattering leads not only to the 
excitation of planetesimal random motions but also to the spatial redistribution of their kinetic 
energy and surface density (see also Paper 11). This is equivalent to the action of the effective 
viscosity in the planetesimal disk which tries to smooth out the gradients of the surface density 
and velocity dispersions S and Sz- As we have mentioned in §2.3 the magnitude of the effect 
produced by the viscosity is typically the same as that of the direct excitation of planetesimal 
velocities by planetesimal-planetesimal scattering. Planetesimal viscosity is very large initially 
when planetesimal velocities are small. As a result the spatial profiles of A^, S and Sz are very 
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Fig. 8. — The same as Figure 7 but for Narn = 10^, p = 10~^. 



smooth (red or blue curves). Later on gravitational effects of the embryo increase as well as the 
planetesimal velocities; this acts to diminish the role of viscosity — spatial distributions start to 
exhibit sharp features driven by the embryo's gravity (green and magenta curves). Finally, deep 
in the asymptotic regime (black curve) where s is very large viscosity is very small: it cannot even 
smooth out strong gradients of S and Sz at the edges of the region occupied by the planetesimals 
on the horseshoe orbits (planetesimal heating inside this zone is very weak). This justifies the 
validity of the concept of the "heated" zone: when embryo dominates nearby disk dynamics, the 
planetesimal viscosity which couples the excited region to the distant parts of the disk becomes 
relatively unimportant, and the heated region can be considered as dynamically isolated from the 
rest of the disk. 

The general picture of the embryo-disk evolution outlined in Figures 3-6 is in excellent agree- 
ment with the predictions of our simple qualitative analysis presented in §2. We have attempted to 
provide some quantitative comparisons by fitting analytical results of 2.2 and 2.3 to the boundaries 
of various regions in Figures 3-5. For example, using the transition between the regions where em- 
bryo and planetesimals separately dominate the disk dynamics (where the color of the dots changes 
in our plots) we can estimate the approximate position of a critical point (s*, A4*) defined by (23): 



« 40 {Narnf^^ ^ 3.4 {Narn)^/^. 



(40) 




In making this estimate we have neglected the variation of logarithmic factor with Narn and 
treated it as a constant which is absorbed in the coefficient in (40). We forced this point to 
lie on the line Ai = (thin solid line) separating shear- and dispersion-dominated regimes of 
the embryo-planetesimal interaction. In the same spirit we have found that the position of the 
separatrix between the regions of planetesimal and embryo dominance over the disk dynamics 
■M.dd{s) [see equation (22)] can be roughly fitted by 

1/2 o 



Mddis) 



M 



5/3 



InAp. 



.1/2 



A. 



1 + 



(41) 



This relation is shown by the long-dashed curve; one can see that within a factor of several it 
correctly predicts the transition from planetesimal to embryo's dominance of the disk dynamics 
although the agreement gets worse as one moves towards the shear-dominated regime [the con- 
tinuation of this separatrix into the shear-dominated region is drawn using (21)]. We have kept 
the logarithmic dependence on s in (41) since it slightly improves agreement with the numeri- 
cal results. Our estimate of M.dd agrees with the corresponding condition suggested by Ohtsuki & 
Tanaka (2002), see their equation (39). They have the same functional dependencies^ and numerical 



^Note that the estimate of Ohtsuki & Tanaka (2002) mistakenly lacks a factor InAe (meaning that embryo- 
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Fig. 10. — Growth of the disk velocity dispersion s in the heated zone as a function of time r 
for Narn = 400, p = 0.004. Sequences of dots display s(t) for different initial conditions. The 
dotted line has a slope of 1/4 and is intended to illustrate the initial growth of s due to the 
planetesimal-planetesimal scattering alone [see (30)]. The dashed line has a slope of 5/9 and shows 
the asymptotic behavior of s. 

coefficients of the same order of magnitude. 

One can see that the evolutionary tracks which were calculated without taking account of the 
embryo's effect on the disk dynamics (thin solid curves) coincide very well with those which include 
embryo's gravity while the disk evolution is dominated by planetesimal-planetesimal scattering. 
This is expected because when M. < Mdd the embryo cannot perturb the disk surface density 
or appreciably affect planetesimal velocities even in its own vicinity. As a result the disk stays 
homogeneous and heats up only by planetesimal encounters — exactly the situation for which 
the curves with a "passive" embryo were computed. However, as soon as the embryo takes over 
the control of the disk dynamics the evolutionary track with an "active" embryo departs from 
that of a "passive" embryo in agreement with equation (33): disk heating speeds up under the 
action of the embryo's gravity while the growth of embryo's mass slows down as a result of the 



planetesimal interaction is in the dispersion-dominated regime) while still retaining In Ae . 
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Fig. 11. — The same as Figure 7 but for Narn = 10^, p = 10~^. 

weaker gravitational focussing of incoming planetesimals and, correspondingly, smaller collision 
cross-section. 

Finally, in the late stages of its evolution the embryo-planetesimal system settles onto an 
attractor. There the dependence of on s can be fitted by a power law with the index 9/5 [see 
equation (28)]. For a given Narn and p all evolutionary tracks independent of their starting Mq 
and So end up on the same curve exactly as envisaged by our analysis in §2.3. Using numerical 
results for different Narn and p this curve can be approximated by [see (28)] 

MM-em{s) ^ 90 [p{NarH)f'''s''/'>. (42) 

This relation is shown by the short-dashed line in Figures 3-5. It is clearly an attractor to which 
all possible initial configurations finally converge. At some point this solution reaches the isolation 
mass Mis which is displayed by the short dash-dotted line — given by (38) in the shear-dominated 
regime and by (39) in the dispersion-dominated regime. It is important to stress that Mis is 
reached in the asymptotic regime of the embryo's evolution and not on the runaway stage. This 
has profound implications for the formation timescale of planetary embryos (see §4). 

We have also checked the agreement between the predictions of §2.3 and our numerical results 
for the location of the separatrix between the regions where the embryo grows faster or slower than 





planetesimals heat up the disk. This is the hne where M. = M.M-pi{s) or Tpi = T^m (see §2.3). 
In numerical calculations we set this transition at the point where the slope of evolutionary tracks 
in s-M coordinates is equal to 1, i.e. where d\n^A/d\ns = 1. The resulting curves for different 
Narn and p are well fitted by [see equation (27)] 

Mm-pi{s) ^ (43) 

where the logarithmic factor InAg is again absorbed into a constant coefficient. 

Note that the values of p and Narjj we have used in our calculations typically do not lead 
to the crossing of Mm-pi{s) with A^M-m(s) in the dispersion-dominated regime [meaning that 
p > po, see equation (29)]. Nevertheless, our numerical results are still in good concord with the 
analytical results of §2.3 which were derived assuming that these curves do cross in the dispersion- 
dominated regime [i.e. p < po]. This is explained by the argument put forward in the end of §2.3: 
values of Narn and p that we have considered are not very extreme {Narn is not large enough and 
p is not small enough) and, as a result, evolutionary tracks do not explore regions of s-^A space 
characteristic for the shear-dominated case. They only dwell for some time in the intermediate 
zones where transitions from one asymptotic regime to the other occur. Of course, if one were to 
explore e.g. larger values of Naru some differences from the case studied in §2.3 would appear. 
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Fig. 13. — Instantaneous surface density of planetesimals which can be accreted by the embryo 
Ninst as a function of time. Curves for several different values of initial embryo's mass are shown 
(labeled on the figure). Note the sharp drop of Ninst after embryo starts to dominate the dynamics 
of the heated zone. 

We expect though that this would not change our principal conclusions. 



3.3. Time dependence of M. and s. 

In Figures 7-9 we plot the dependence of embryo's mass on time for the evolutionary tracks 
shown in Figures 3-5. Here one can see better that disk evolution with the "passive" embryo (thin 
sohd curves) really leads to a runaway — mass of the embryo grows to extremely large values in a 
finite time. We have found that mass growth on this runaway stage can be approximately described 
by equation (30) with r^un roughly (within a factor of 2) given by 

0-0^ (44) 



Evolutionary paths corresponding to the more realistic case of an "active" embryo depart from 
the runaway tracks when the embryo starts dominating disk dynamics in its vicinity, just as in 
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Figures 3-5. They smoothly switch to a behavior where is a power law of time [similar behavior 
was obtained in Weidenschilling et al. (1997)]. Although this is consistent with the results of our 
qualitative analysis the exact value the power law index does not coincide with the value of 9/7 
predicted by equation (34). In fact it is closer to 1 so that the embryo's mass grows linearly with 
time. We have found that the time dependence of mass can be roughly described by the following 
relation: 

M{t) ^ 1300 (7Var//)^-^V V. (45) 

There is a simple explanation of this apparent discrepancy: when the embryo dominates the 
disk heating it induces strong nonuniformities in the distribution of the planetesimal surface density 
around it. As a result the instantaneous surface density of planctesimals at the embryo's location 
which determines its accretion rate decreases with time. We illustrate this point in Figure 13 where 
we plot the instantaneous surface density of planetesimals near the embryo Ninst normalized by its 
initial value as a function of time r for Narn = 10^, p = 10~^ and several different values of initial 
mass A^o- Values of Ninst are computed from spatial profiles of N and S using the corresponding 
conversion found in Paper II; in calculating Ninst we exclude the contribution of the horseshoe 
orbits because they cannot approach the embryo very closely and be accreted. One can see that 
Ninst begins to vary strongly only after the embryo starts dominating the disk dynamics (where 
the color of dots changes to red); in this regime Ninst can drop by a factor of 10 or more because 
of the spatial redistribution of planetesimals in the disk. This is equivalent to a decrease of Narn 
which would otherwise be constant. As a result, the accretion rate diminishes and embryo's mass 
grows slower than equation (34) would predict. The small variations of Ninst that occur in the 
stage when planetesimal dynamics is determined only by planetesimal scattering (blue dots) are 
explained by our exclusion of planetesimals on horseshoe orbits^: the spatial extent of the horseshoe 
region depends on the eccentricity and inclination of the planetesimals (see Paper III for details) 
which vary with time; this causes Ninst to change. 

Our simple theory presented in §2.3 does not account for the local decrease of the planetesimal 
surface density as M the planetesimals are accreted onto the embryo. This is a considerable 
drawback which requires a more sophisticated treatment. We will not pursue this subject in this 
paper, merely note on its importance. Luckily it docs not affect the general picture of the evolution 
of embryo-disk system described in §2.3, up to the isolation mass Mis- 
Note a very strong increase of the embryo's growth timescale in Figures 7-9 as one moves 
towards the outer edge of the nebula. The dimcnsionless timescale (in units of r) grows because it 
strongly depends on parameter p characterizing embryo's accretion cross-section and p diminishes 
when the distance from the central star a increases. The timescale measured in physical units is 
additionally lengthened by the growth of the planetesimal synodic period tgyn with a: tgyn oc a^/^. 



*In this regime the planetesimal disk is almost homogeneous; if we did not exclude planetesimals on the horseshoe 
orbits Ninat would be exactly constant. 
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As a result the most favorable conditions for the embryo's growth exist in the inner parts of the 
nebula. We will return to this issue in §4. 

In Figures 10-12 we plot the dependence of s on time in the same way as we did it for Ai. 
Initially (when planetesimal encounters control disk dynamics) one would expect s to behave in 
accordance with equation (30) which is exhibited by dotted line. Deviations from this power-law 
solution are caused by additional logarithmic dependence present in equation (30). At some point 
evolution switches to the embryo-dominated regime. We have used Figures 10-12 to determine the 
locations of the corresponding transition points for each track by approximately picking the point 
at which there is a break in the s(r) dependence. As wc have discussed above, this choice of the 
transition points agrees well with other characteristic indicators e.g. splitting of A1(r) for "passive" 
and "active" embryo tracks in 7-9. When r becomes very large system evolves in the asymptotic 
regime and s grows as a power-law of r with index close to 5/9. This particular dependence and 
the fact that A^(t) oc r in the asymptotic regime conspire to produce the dependence of A4{s) 
which is very close to that predicted analytically [see (28) and (42)]. The fact that the power-law 
index of s(t) is different from 5/7 predicted by (34) is again explained by the drop of Ninst in the 
asymptotic regime: since embryo grows slower that expected it heats up planetesimal population 
nearby less efficiently than our analysis of §2.3 predicts. Nevertheless, by the end of our calculations 
planetesimal velocity dispersion in the heated zone exceeds the random velocity in the distant part 
of the nebula typically by a factor of several. As expected, the gravitational influence of the 
embryo considerably speeds up the dynamical evolution of the planetesimals in the heated zone. 
This behavior is in general agreement with the results of multi-zone simulations by Weidenschilling 
et al. (1997). 



4. Discussion. 

4.1. Results and their applications. 

The most important result of our investigation is the elucidation of the role played by the 
coupling between the nonuniform disk dynamics and the embryo's mass growth. We have found 
in agreement with the previous investigations (Ida & Makino 1993; Tanaka & Ida 1997; Paper I) 
that the dynamical excitation of the planetesimal disk by the gravity of growing protoplanetary 
embryo dramatically changes the simple runaway picture of planetary growth, which is derived for 
homogeneous planetesimal disks. Using our results we can also study quantitatively the question 
of how long is required for the embryo to form? 

As we have pointed out in §3, embryo accretes all the material within its feeding zone and 
reaches the isolation only after it settles into the asymptotic evolutionary regime (and not on 
the runaway stage). The time that the embryo spends on this asymptotic stage dominates the 
final timescale of its evolution. Indeed, for all the values of Narn and p explored in §3 the 
runaway timescale is usually shorter by a factor of several than the time needed for the embryo 
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to reach isolation in the asymptotic stage. This might indicate that the estimates of the embryo 
formation timescale produced by the conventional coagulation simulations neglecting the effects of 
inhomogeneous disk evolution could be off by a substantial factor. 

For example, in the case Narn = 400, p = 0.004 typical for the planet formation in the 
terrestrial planet region of MMSN, the time needed for "passive" runaway embryo (thin solid 
curve) with Mg/m = 10 to reach the isolation mass 3 x 10^^ g is 4 x 10^ yr; growth to the 
same mass for an "active" embryo that affects the disk would take ~ 3 x 10^ yr (see Figure 7) 
— almost a factor of 10 increase! This difference is even more striking for larger initial M^/m for 
which runaway timescales arc very short. 

The very long time needed for the embryo's growth which arises in our calculations exacerbates 
the timescale problem for the growth of the rocky cores of the giant planets, because this process 
usually takes at least as long as the gaseous nebula is supposed to last: at 3.6 AU in the MMSN (a 
situation described by p = 10~^ and Narn = 10^ in our computations) the growth of Mg to 2 x 10^''' 
g [the isolation mass given by (38)] takes 10^ yr, while reaching 9 x 10^^ g [isolation mass given 
by (39) which accounts for the expansion of the feeding zone due to the excitation of planetesimal 
random motions] would take 5 x 10^ yr. Not only is this timescale is too long, it is also not 
clear that the isolation mass is large enough to trigger the core instability which would allow giant 
planets to accrete their huge gaseous mass (Mizuno 1980; Bodenheimer & Pollack 1986). Post- 
isolation chaotic dynamical evolution of the embryos (Chambers 2001) is the most likely solution 
of the latter problem, but it would require even more time for building up massive gas-accreting 
planetary cores. Note that the existence of this post-runaway chaotic evolutionary stage would not 
be a problem for the formation of terrestrial planets because (1) timescales are shorter in the inner 
parts of the nebula (see Figure 7) , and (2) terrestrial planets did not accrete large amounts of gas 
which might indicate that they have formed after the gaseous nebula dispersed (which eliminates 
the time constraint related to the nebular dispersal). 

It is true that in our calculations the isolation mass is usually larger than Aiig reached by 
the embryos growing in the runaway regime in homogeneous disks: in the runaway stage Mis is 
given by (38) since it is typically reached in the shear-dominated regime (see Figures 3-5) while for 

the "active" embryos Mis is increased as a result of the strong excitation of planetesimal random 
velocities by the embryo, see equation (39). This simplifies the formation of massive rocky cores 
(at 3.6 AU Mis is increased by almost a factor of 5, see above). Unfortunately, the time required 
to reach larger Mis is longer as well. 

The timescale issue becomes even more challenging for the growth of the outer giant planets 
if they formed at their present locations^: typical formation timescales are ~ 10^ yr (see Figure 9) 
which is much longer than the timescale of the nebular dissipation (Hollenbach et al. 2000). One 



®It was suggested by Thommes et al. (2002) that Uranus and Neptune were formed between the orbits of Jupiter 
and Saturn and were then scattered into their present orbits by these giant planets. 
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possibility to alleviate this problem would be to increase the surface density of planetesimals; this 
is a reasonable assumption beyond the so-called "snow line" (Hayashi 1981; Sasselov &; Lecar 2002) 
where the condensation of ice takes place (at 1 — 3 AU) . 

The fact that the embryo's growth does not stop when the embryo-planetesimal interaction 
starts to dominate disk dynamics and opens a gap is quite interesting. Previous estimates of the 
gap opening mass (Paper I; Ohtsuki &; Tanaka 2002) only considered an embryo with a fixed mass 
and allowed plenty of time for the disk evolution. Under these conditions a gap would form in 
the planctcsimal disk around the embryo when Tgm ~ Tpi and Ai ~ Mddis); the instantaneous 
accretion rate of the embryo (if it were able to increase its mass) would be strongly affected even at 
this rather small value of Ai. This would strongly slow down embryo's accretion at Ai ^ -Mddis)- 
This reasoning essentially assumes that the disk instantaneously adjusts to any changes of the 
embryo's mass. 

The picture turns out to be quite different in reality because the disk cannot quickly adapt to 
the embryo's evolution — on the contrary it is the embryo that initially evolves much faster than 
the disk when Ai > Aidd- A-S a result, during the time needed for the embryo's gravity to clear a 
gap in the disk its mass grows substantially. Although accretion in this part of s-Ai space is not in 
the runaway regime it is still quite fast initially; only after the embryo gains a considerable amount 
of mass the accretion smoothly changes its character to an orderly growth with Tm ~ in the 
region between the curves AA. = Aidd and Ai = AiM-em- Thus, accretion sets into an orderly 
mode not when Tg^ Tpi (as analysis with a fixed embryo mass would assume) but only when 
Tf,jn Tm which implies larger critical mass of the embryo-'^'^ . 

Strong increase of the planetesimal velocity in the heated zone compared to that at infinity can 
be an important issue for the planetesimal growth in this region. Indeed, planetesimal coagulation 
in this part of the disk would almost certainly proceed in the orderly regime (in contrast to the 
embryo's accretion): planetesimal velocities should be so high that the gravitational focussing is 
negligible and accretion cross-section is independent of planetesimal velocities [see equation (3)]. 
Planetesimal agglomeration in the excited zone would then be so slow that the growth of more mas- 
sive bodies in the heated region might essentially be chocked (planetesimal fragmentation would 
only strengthen this conclusion). New planetary embryos would be able to emerge only far enough 
from the preexisting ones where the planetesimal dynamics is not strongly affected by the gravita- 
tional perturbations of massive bodies. 

It is also interesting to note that gravitational focussing is always important for the embryo's 
accretion (unlike the case of orderly growth, see §1). Focussing is weak when the planetesimal 
velocities are larger than the escape speed from the embryo's surface or when S > 3> 1 

or s > 

MV3p-i/2 [ggg (A21)]. In Fi gures 3-5 we plot this restriction with a long-dash-dotted 
line (5 = p-V2). One can see that evolutionary tracks never penetrate into this region: although 



As we have mentioned above this increase is not big enough to resolve the timescale issue. 
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initially S increases towards p^^^^ as planetesimal velocity grows, the rapid increase of the embryo's 
mass caused by the runaway accretion makes S drop significantly (sometimes bringing embryo- 
planetesimal interaction into the shear-dominated regime) when the condition Tm < Tpi is fulfilled, 
see §2.3 and equation (43). This implies that the embryo's accretion cross-section is always strongly 
increased over its geometric value by the gravitational focussing-*--^ . An important conclusion which 
can be drawn from this is that the growth of the embryo never proceeds in the orderly fashion 
with weak gravitational focussing, as proposed by Safronov (1972). Instead, it starts as a runaway 
evolution (Wetherill & Stewart 1989) which then switches to an oligarchic growth similar to that 
predicted by Kokubo & Ida (1998). 



4.2. Limitations of the analysis. 

The importance of our conclusions relies on the applicability of our simple model to the realistic 
protoplanetary disks. There is a lot of complications which can potentially affect our results. The 
assumption of a single mass population is certainly one of them. Our whole analysis relies on the 
possibility to compress the details of the planetesimal mass spectrum into properties of a single 
population with a unique characteristic mass. This is only possible if the estimates of such averaged 
mass characterizing separately the embryo's accretion and disk dynamical evolution agree with each 
other. If the distribution of masses and random velocities is such that this condition is not fulfilled 
(i.e. if embryo's accretion and disk dynamics are dominated by different parts of the planetesimal 
mass spectrum) then our consideration would fail and one would need to resort to a more intricate 
analysis. 

Inclusion of a spectrum of planetesimal masses would also allow dynamical friction to redis- 
tribute the random energy of epicyclic motion between planetesimal populations with different 
masses (Stewart &; Wetherill 1988). This might lead to an interesting situation when the embryo's 
interaction with different populations proceeds in different dynamical regimes — some will be in the 
shear-dominated regime while some in the dispersion-dominated. This can introduce some quanti- 
tative corrections to our results but we would expect the general picture of the planetary growth 
outlined here to remain unchanged: we have demonstrated in §2.3 that even with a single mass 
population of planetesimals we can still adequately reproduce the major features of the embryo's 
runaway growth. 

Another potential worry is the evolution of the mass spectrum. This leads to the evolution 
of the average planetesimal mass m and affects Narn as a result. Note that changes of m can 
be caused not only by planetesimal coagulation which always tends to increase m but also by the 
change of the shape of the planetesimal mass spectrum driven by the embryo's presence. The last 
possibility might in principle lead to a decrease of m if the control over the embryo-disk evolution 



This also means that embryo's erosion by planetesimals colliding with it is unimportant. 
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switches to a lower-mass part of the planetesimal spectrum. Evolution of m is probably not very 
important initially because in the runaway regime the embryo grows faster than planetesimals, see 
equation (4). But later on when embryo dominates the disk dynamics and the system settles into 
the asymptotic regime with very slow accretion the planetesimal mass evolution might lead to some 
new effects. In addition, the growth of m is tightly coupled to the damping of the planetesimal 
epicyclic motion by inelastic scattering because both damping and coagulation are consequences of 
the same process — direct planetesimal collisions. Preliminary calculations demonstrate however 
that the decrease of planetesimal velocities by inelastic collisions would not allow embryo's accretion 
to continue in the runaway regime until the isolation mass is reached — at some point it would 
still switch to an orderly regime where M grows as some power law of time, although this orderly 
growth will be somewhat faster than that predicted by our analysis of §2.3 [see (34)]. Planetesimal 
fragmentation can also become important when relative planetesimal velocities exceed the escape 
speed from the planetesimal surface (roughly when s > One can extend our analysis and 

numerical calculations to include at least the most important of these effects into account and we 
are going to do this in the future. 

Our consideration has also been restricted to studying only a single embryo in the disk. In 
reality there would probably be many of them, at least initially. Since their masses grow by several 
orders of magnitude the Hill radii of neighboring embryos (as well as their excited regions) would 
overlap at some point. This would homogenize the disk to some extent but it would also slow 
down the accretion even more because disk would be heated by several embryos at a time. Finally 
gravitational perturbations between embryos would lead to their collisions and agglomeration. This 
question certainly deserves more attention both from the theoretical and from the numerical points 
of view. One should also be concerned with the possibility of embryo's migration which can in 
principle replenish the feeding zone of the embryo (Tanaka & Ida 1999; Ida et al. 2000). 

Finally, gas drag can be important for the planetesimal dynamics especially for small plan- 
etesimal masses (Adachi et al. 1976). Its effect is similar to that of inelastic collisions — it damps 
epicyclic motions of planetesimals and leads to somewhat faster accretion by the embryo in the 
asymptotic regime. Gas drag can naturally be incorporated into our analytical apparatus of Paper 
II and its effects will be explored later on. 

5. Conclusions. 

In this paper we have self-consistently studied the dynamical evolution of a planetesimal disk 
coupled to the growth of a single massive protoplanetary embryo. We are able to demonstrate that 
the evolution of the embryo-disk system proceeds in two consecutive stages: 

• It starts with a rapid runaway growth of the embryo, during which its mass is not large enough 
to affect the disk dynamics and the spatial distribution of planetesimals. On a rather short 
timescale (comparable to that arising in conventional coagulation simulations) the embryo 
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reaches the mass at which it takes over the control of the disk heating and runaway growth 
stops. 

• After a short transient stage, the embryo-disk system settles into the asymptotic regime in 
which the timescale of the embryo's mass growth is comparable to the timescale on which the 
disk is heated by the embryo's gravity. 



During the last stage, the embryo dynamically excites planetesimal epicyclic motions in its vicinity 

and repels the planetesimals forming a depression in the surface density of the guiding centers (see 
Figure 6). This effect produces a negative feedback for the embryo's accretion rate for two reasons: 
rapid growth of the random velocities of planetesimals leads to weaker gravitational focussing while 
the decrease of the instantaneous surface density (see Figure 13) reduces the amount of material 
available for accretion. As a result, growth of the protoplanetary embryo proceeds slower than in 
homogeneous planetesimal disks. This implies that conventional "particle-in-a-box" coagulation 
simulations (Wetherill & Stewart 1989; Kenyon & Luu 1998; Inaba et al. 2001) might be underes- 
timating the timescale of the protoplanetary embryo's growth, sometimes by rather large factors 
(up to ~ 10). Multi-zone coagulation simulations (Spaute et al. 1991; Weidenschilling et al. 1997) 
should give a more reliable description of the protoplanetary disk evolution. 

We have presented OTir equations in a dinicnsionlcss form which has allowed us to uncover a very 
important property of the problem we consider: its outcome depends on only two dimensionless 
parameters — the number of planetesimals inside the annulus of the disk with width equal to 
the planetesimal Hill radius rn, Narn, and the ratio of the physical to Hill radii of a solid body 
(planetesimal or embryo) p. All astrophysically relevant characteristics of the system (m, a, Mc, S, 
p) are combined into these two parameters which greatly simplifies the exploration of the parameter 
space. We have numerically studied in §3 the evolution of the system for a number of pairs of Narn 
and p; combined with our analytical developments in §2 these results allow us to formulate a set 
of simple scaling relations (40)- (45) which can be used to predict the embryo-disk evolution for 
different sets of Narn and p. 

To summarize we can say that the embryo's growth starts in the runaway regime, but then 
switches into the "oligarchic" growth of Kokubo k, Ida (1998). We have demonstrated that the 

embryo's accretion never proceeds in the orderly regime with weak gravitational focussing suggested 
by Safronov (1972) — planetesimal velocities are always smaller than the escape velocity from the 
embryo's surface. The embryo formation timescale and some details of the embryo-disk evolution 
are subject to a number of uncertainties related to the simplicity of our model, as we discussed in 
§4. Nevertheless, we believe that our major conclusions are robust: the results of our qualitative 
analysis of §2 and numerical calculations presented in §3 suggest that the dynamical interaction 
between the protoplanetary embryo and planetesimal disk is a very important issue which should 
certainly be addressed when realistic coagulation simulations are performed. The general agreement 
between the results of our simple analysis of §2 with more accurate numerical calculations of §3 
encourages the studies of more complex and more realistic problems in the same spirit. In the future 
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we are going to improve our treatment of the coupled planetesimal disk and embryo evolution by 
relaxing simplifying assumptions employed in this paper. We are also going to include additional 
mechanisms which are important in modelling of planetary formation — effects of multiple embryos, 
migration, gas drag, etc. 

I am grateful to Scott Tremaine for his advice, and for reading this manuscript and making 
a lot of useful suggestions. The financial support of this work by the Charlotte Elizabeth Procter 
Fellowship and NASA grant NAG5-10456 is thankfully acknowledged. 
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A. More accurate equations of evolution. 

Here we describe the system of equations we are going to use to follow the evolution of plan- 
etesimal surface density, eccentricity and inclinations. We derive these equations for the general 
case of a distribution of planctesimal masses for future applications. For this reason here N is a 
function of planctesimal mass"*^^ m. 

We will use the following notation. First of all, we introduce the dimensionless radial distance 
from the embryo's orbit H, normalized by its Hill radius: 

H = Aa/RH, RH = aixl'^, iJLe = MjM^, (Al) 

where Aa is the dimensional radial distance from the embryo's orbit (equivalent to ha in the 

1/3 

notation of Paper H). We also scale all eccentricities and inclinations by Rh/cl = A*e ■ This gives 
us new variables S and Sz defined by equation (35). When doing this we have to keep in mind that 
S and Sz can vary not only because eccentricity and inclinations (t^ and cTj vary with time but also 
because /Xe can change with time (as a consequence of planetary mass growth). This variation can 
be easily taken into account by noticing that 

We will characterize mass distribution of planetesimals by some fiducial mass rrii, and introduce 
a dimensionless mass z such that 

z = m/m^ = n/fi^, II = m/Mc, /x* = rUi^/Mc. (A3) 

Then M = Mg/m* Associated with this mass is a fiducial surface number density A?* = Tto/rrii,, 
where Sq is the total surface mass density of the disk. Following Papers H &: HI we assume both 
N{Tn) and Sq to be scaled by a^ (so that N{m) = N{z) is dimensionless, unlike §2). Finally we 
define a new surface number density function 

oo 

g{z, H) = = ml^, J g{z, H)zdz = 1. (A4) 



We will use the notation gi{H) for the surface number density of planetesimals of mass rrii (or 

In transforming the planetesimal-planetesimal part of the evolution equations we will take into 
account that scattering coefficients in the dispersion-dominated regime have a rather specific form 
given by equations (82)- (93) of Paper II; we need to remember that expressions for these coefficients 



Unlike §2 where we have studied only a single mass population. We go back to this simpler case later on. 
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are functions of C7e,i normalized by Hill radius of planetesimal-planetesimal interaction rn, which is 
different from S and S^- With this in mind, using definitions (A1)-(A4) and equations (49)-(55) of 
Paper II and (13)-(15), (24)-(26), (B2) of Paper III we can finally represent equation of evolution 
of quantity F = {g,gs^ ,gs1} in the following general form: 



pi 



dF 



hot 



^0 



cold- 



2 dlnM 
~ 3^ dr 



(A5) 
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with a new dimensionless time r = (3/47r)Ot/// = t/tgyn [see (14)]. The interpolating functions 
6i = Q{S, Sz) and 62 = I — &{S, Sz) provide smooth matching between the shear- and dispersion- 
dominated regimes of embryo-planetesimal interaction (marked with subscripts "cold" and "hot" , 
see Paper III). Function @{S,Sz) has the following properties: @{S,Sz) 1 as S,Sz ^ 00, and 
0(5, 5^) ^ as 5, 5*2 — 0. The last term in (A5) is present only in equations for eccentricity and 
inclination evolution. 



For the surface density of planetesimals of type 1, F = gi and we have 
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where 



T^{H) = 2 j dz2Z2{zi + Z2) j dHi\Hi\{Ah)g2{H - Hi), 

-00 

00 00 

T^{H) = j dz2zl j dHi\Hi\{{/\hf)g2{H-Hi), 
-00 

where Z2 is used as a variable of integration over the planetesimal mass spectrum. 
For eccentricity evolution, F = 2giSi and 
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where 



oo 

+ J P{Ho,H)gi{Ho)(2sliHo) + ^A{esc)\Ho)^\Ho\dHo, (A13) 
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Tl{H) = 2 j dz2Z2{zi + Z2) j dHilHMH - Hi) 
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Si = Si{H), S2 = S2{H-Hi). 
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Different terms in the equation of evolution of 2giSli look analogous to (A11)-(A12) and 
(A14)-(A16) if we replace 5 with Sz but the shear-dominated part reads (see Paper III) 
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2gi{H)Sli{H)\H\-2 / P{Ho,H)gi{Ho)Sli{Ho)\Ho\dHo. (A18) 



In (A8), (A13), (A18) Hq is an integration variable having the meaning of initial difference in 
semimajor axes of interacting bodies and 



P{Ho, H) = S H-Ho- AhsciHo) 



(A19) 



where hsdHo) is some function for which the analytical prescription based on the results of numer- 
ical calculations exists (see Petit & Henon 1987). Scattering coefficients {Ah), {{Ah)'^), {e^Ah), 
((e • Ae)Ah), {e'^{Ah)'^), ((Ae)2), (e • Ae) in these equations are functions of 5, Sz, e.g. 

(A^) = 4^e--V(25.^)lt;0 (H/S^H/S. 
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H ^\2V2 2V2 



A = A^i/3^^(c,52 + ^^^2^) 1^ ^ ^52 + 52, Szr = \lsli + Sl^, (A20) 
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where ci and C2 are some constants (their numerical values are fixed approximately in Paper III), 
and function Up is defined by equation (95) of Paper II. Scattering coefficients describing embryo- 
planetesimal interaction are analogous to (A20) but have A = (ciS^ + C2S'^^)Szr/{zi + 2:2)- 

For the growth of the planetary embryo's mass we can write using equations (42)- (44) of Paper 

III: 

00 





00 



(A21) 



where functions vri = 11(5', Sz) and = 1 — 11(5", Sz) are used to interpolate between the shear- and 
dispersion-dominated regimes of embryo-planetesimal interaction. Function 11(5', Sz) has properties 
analogous to the properties of @{S,Sz) (see Paper III). The quantity 5*"** in the second term in 
brackets is the instantaneous surface density of planetesimals at the embryo's location, which can 
be calculated using the conversion between N and AT'"** described in Paper III. The dimensionless 
parameter p was defined in (14). The first terms in the square bracket describe the mass accretion 
rate in the dispersion-dominated and shear-dominated regime, respectively. 

When a single mass population of planetesimals is considered one should use the planetesimal 
mass m as a fiducial mass and set g{z, H) = 5{z ~ l)f{H), where f{H) is a function of distance 
H only. Then the integration over Z2 in (A9), (AlO), (A14)-(A16) trivially goes away. In this work 
we have only considered a single mass population of planetesimals. 

One can easily notice that in agreement with the results of §2 there are only two free parameters 
present in equations (A5)-(A21): p — ratio of physical radius to the Hill radius of any object, and 
Ni,ii\/^ — measure of the disk surface number density equivalent to Narn of §2 (in §2 AT is 
dimensional). When the distribution of planetesimal masses exists the shape of the mass spectrum 
is also a prescribed input function (which can vary in time as a result of coagulation). Equations 
(A5)-(A21) constitute a closed set of equations describing the behavior of the embryo-disk system. 
Their numerical solutions are presented in §3. 



